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<^ ! Abstract 



Magneto-Acousto-Electric Tomography (MAET), also known as the Lorentz force 
or Hall effect tomography, is a novel hybrid modality designed to be a high-resolution 
alternative to the unstable Electrical Impedance Tomography. In the present paper we 
analyze existing mathematical models of this method, and propose a general procedure 
_C ■ for solving the inverse problem associated with MAET. It consists in applying to the 

data one of the algorithms of Thermo- Acoustic tomography, followed by solving the 
Neumann problem for the Laplace equation and the Poisson equation. 

For the particular case when the region of interest is a cube, we present an explicit 
series solution resulting in a fast reconstruction algorithm. As we show, both analyt- 
ically and numerically, MAET is a stable technique yilelding high-resolution images 
■ even in the presence of significant noise in the data. 

CO 

q . 

od : Introduction 

o . 

Magneto-Acousto-Electric Tomography (MAET) is based on the measurements of the elec- 
trical potential arising when an acoustic wave propagates through conductive medium placed 
in a constant magnetic field [26,34]. The interaction of the mechanical motion of the free 
charges (ions and/or electrons) with the magnetic field results in the Lorentz force that 
pushes charges of different signs in opposite directions, thus generating Lorentz currents 
within the tissue. The goal of this technique coincides with that of the Electrical Impedance 
Tomography (EIT): to reconstruct the conductivity of the tissue from the values of the 
electric potential measured on the boundary of the object. EIT is a fast, inexpensive, and 
harmless modality, which is potentially very valuable due to the large contrast in the con- 
ductivity between healthy and cancerous tissues [5,6,9]. Unfortunately, the reconstruction 
problems arising in EIT are known to be exponentially unstable. 

MAET is one of the several recently introduced hybrid imaging techniques designed 
to stabilize the reconstruction of electrical properties of the tissues by coupling together 
ultrasound waves with other physical phenomena. Perhaps the best known examples of 
hybrid methods are the Thermo-Acoustic Tomography (TAT) [16] and the closely related 
Photo-Acoustic modality, PAT [15, 29]). In the latter methods the amount of electromagnetic 
energy absorbed by the medium is reconstructed from the measurements (on the surface 
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of the object) of acoustic waves caused by the thermoacoustic expansion (see e.g. [17, 
33]). Another hybrid technique, designed to overcome shortcomings of EIT and yield stable 
reconstruction of the conductivity is Acousto- Electric Impedance Tomography (AEIT) [36]. 
It couples together acoustic waves and electrical currents, through the electroacoustic effect 
(see [25]). Although AEIT has been shown, both theoretically and in numerical simulations, 
to be stable and capable of yielding high-resolution images [3,8,18,19], the feasibility of 
practical reconstructions is still in question due to the extreme weakness of the acousto- 
electric effect. 

In the present paper we analyze MAET which also aims to reconstruct the conductivity 
in a stable fashion. In MAET this goal is achieved by combining magnetic field, acoustic 
excitation and electric measurements, coupled through the Lorentz force. The physical 
foundations of MAET were established in [26] and [34]. In particular, it was shown in [26] 
that if the tissue with conductivity cr(x) moves with velocity V(x,t) within the constant 
magnetic field B, the arising Lorentz force will generate Lorentz currents Jl(x, t) whose 
intensity and direction are given (approximately) by the following formula 



Originally [26,34] it was proposed to utilize a focused propagating acoustic pulse to 
induce electric response from different parts of the object. In [13] wavepackets of a certain 
frequency were used in a physical experiment to reconstruct the current density in a thin 
slab of a tissue. Similarly, in [4] the use of a perfectly focused acoustic beam was assumed 
in a theoretical study and in numerical simulations. However, in the above-quoted works 
accurate mathematical model (s) of such beams were not presented. Moreover, the feasibility 
of focusing a fixed frequency acoustic beam at an arbitrary point inside the body in a fully 
3D problem is problematic. In a theoretical study [30] the use of plane waves of varying 
frequencies was proposed instead of the beams. This is a more realistic approach; however, 
the analysis in that work relies on several crude approximations (the conductivity is assumed 
to be close to 1, and the electric field is approximated by the first non-zero term in the 
multipole expansion). 

To summarize, the existing mathematical models of measurements in MAET are of ap- 
proximate nature; moreover, some of them contradict to others. For example, it was found 
in [13] that if one uses a pair of electrodes to measure the voltage (difference of the potentials) 
at two points a and b on the boundary of the body, the result is the integral of the mixed 
product of three vectors: velocity V, magnetic induction B and the so-called lead current 
J a b (the current that would flow in the body if the difference of potentials were applied 
at points a and b). The approximate model in [30] implicitly agrees with this conclusion. 
However, in [4] it is assumed that if the pulse is focused at the point x, the measurements 
will be proportional to the product of the electric potential u(x) and conductivity a(x) at 
that point. This assumption contradicts the previous models; it also seems to be unrealistic 
since potential u(x) is only defined up to an arbitrary constant, while the measurements are 
completely determined by the physics of the problem. 

In the present paper we first derive, starting from equation (1), a rigorous and sufficiently 
general model of the MAET measurements. Next, we show that if a sufficient amount of data 
is measured, one can reconstruct, almost explicitly and in a stable fashion the conductivity 
of the tissue. For general domains the reconstruction can be reduced to the solution of the 
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inverse problem of TAT followed by the solution of the Neumann problem for the Laplace 
equation, and a Poisson equation. In the simpler case of a rectangular domain the recon- 
struction formulae can be made completely explicit, and the solution is obtained by summing 
several Fourier series. In the latter case the algorithm is fast, i.e. it runs in 0(n 3 logn) float- 
ing point operations onanxnxn grid. The results of our numerical simulations show that 
one can stably recover high resolution images of the conductivity of the tissues from MAET 
measurements even in the presence of a significant noise in the data. 



1 Formulation of the problem 

Suppose that the object of interest whose conductivity cr(x) we would like to recover is 
supported within an open and bounded region Q with the boundary dfl. For simplicity we 
will assume that a(x) is smooth in Q, does not approach 0, and equals 1 in the vicinity of 
dfl; the support of a(x) — 1 lies in some C VL and the distance between and dVt is 
non-zero. The object is placed in the magnetic field with a constant magnetic induction B, 
and an acoustic wave generated by a source lying outside Q propagates through the object 
with the velocity V(x,t). Then the Lorentz force will induce Lorentz currents in Q given by 
equation (1). Throughout the text we assume that the electrical interactions occur on much 
faster time scale than the mechanical ones, and so all currents and electric potentials depend 
on t only as a parameter. In addition to Lorentz currents, the arising electrical potential 
u(x,t) will generate secondary, Ohmic currents with intensities given by Ohm's law 

Jo{x,t) = a(x) J \7u(x,t). 

Since there are no sinks or sources of electric charges within the tissues, the total current 
Jl(x, t) + Jo(x,t) is divergence- free 

V • (J L + Jo) = 0. 

Thus 

V • (tV« = -V • (aB x V) . (2) 

Since there are no currents through the boundary, the normal component of the total current 
Jl(x, t) + Jo(x,t) vanishes: 

^-u(z) = -(B x Viz)) ■ n(z), z e dQ, (3) 
on 

where n(z) is the exterior normal to dfl at point z. 

We will assume that the boundary values of the potential u(z,t) can be measured at all 
points z lying on dQ. More precisely, we will model the measurements by integrating the 
boundary values with a weight I(z) and thus forming measuring functional M defined by 
the formula 

M(t) = J I(z)u(z,t)dA(z), (4) 
dn 



3 



where A(z) is the standard area element. Weight I(z) can be a function or a distribution, 
subject to the restriction 

J I(z)dA(z) = 0. 
an 

In particular, if one chooses to use I(z) = 8{z — a) — 8(z — b), where S(-) is the 2D Dirac 
delta-function, then M models the two-point measuring scheme utilized in [13] and [30]. 

In order to understand what kind of information is encoded in the values of M(t) let us 
consider solution wj(x) of the following divergence equation 

V • ffVw/(x) = 0, (5) 

■^-wr(z) = I (z) , zedtt. (6) 
on 

(To ensure the uniqueness of the solution of the above boundary value problem we will 
require that the integral of wj{z) over Vt vanishes.) Then wj(x) equals the electric potential 
that would be induced in the tissues by injecting currents / (z) at the boundary dQ. Let us 
denote the corresponding currents <jVwi{x) by Ji{x): 

Ji{x) = aVwi{x). (7) 

Let us now apply the second Green's identity to functions u(x, t), wj(x), and a(x): 



J [wjV ■ (<tVu) - wV • {aVwj)]dx = J a 



d d 

Wj—U - U—Wj 

on on 



dA(z). (8) 



By taking into account (3), (2), (5), and (6), equation (8) can be simplified to 
- f wjV ■ (aB xV)dx= [ awj^-udA(z) - M(t). 



dn 



n on 



Further, by integrating the left hand side of the last equation by parts, and by replacing -j^u 
with expression (3) we obtain 

J aVwj .(BxV)dx-J wja(B x V) ■ ndA(z) = - J awj(B x V) ■ ndA(z) + M(t) 
n an an 

or 

M(t) = - J aVwj -(BxV)dx = - J J/ • (B x V) dx = B ■ J J^x) x V(x, t)dx (9) 
n n n 

This equation generalizes equation (1) in [13] obtained for the particular case I(z) = S(z — 
a)-S(z-b). 

It is clear from equation (9) that MAET measurements recover some information about 
currents Ji{x). In order to gain further insight let us assume that the acoustical properties 
of the medium, such as speed of sound c and density p, are approximately constant within 
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Q. (Such approximation usually holds in breast imaging which is one of the most important 
potential applications of this and similar modalities). Then the acoustic pressure p(x,t) 
within Q satisfies the wave equation 

1 d 2 

Additionally, p(x,t) is the time derivative of the velocity potential <f(x,t) (see, for exam- 
ple [10]), so that 



V(x,t) = -V<p(x,t), 
P 

P(x,t) = —<p(x,t). 

Velocity potential (p(x,t) also satisfies the wave equation 

1 d 2 



c 2 dt 2 



(p(x, t) = A(p(x, t). 



Now, by taking into account (10), equation (9) can be re-written as 

M(t) = -B-[ JAx) x Vip(x,t)dx 
P J 

Further, by noticing that 

V x (tpJi) = —Ji x V(p + tpV x J/ 



we obtain 



M(t) = -B ■ 



= -B ■ 

P 



- / V x (<fj 1 )dx+ / (p(x,t)V x J I (x)dx 



(p(z,t)Ji(z) x n(z)dA(z) + J ip(x,t)V x J I (x)dx 



(10) 



(11) 



In some situations the above equations can be further simplified. For example, if at some 
moment of time t velocity potential (p(x,t) vanishes on the boundary dfl, then the surface 
integral in (11) also vanishes: 



M(t) = ~ p B - J V?(M)V x J T (x)dx. 



(12) 



Similarly, if boundary dQ is located far away from the support of inhomogeneity of cr(x), the 
surface integral in (11) can be neglected, and we again obtain equation (12). 
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Equation (11) is our mathematical model of the MAET measurements. Our goal is 
to reconstruct from measurements M(t) conductivity cr(x) — by varying, if necessary B, 
ip(x, t), and I(z). 

Our strategy for solving this problem is outlined in the following sections. However, some 
conclusions can be reached just by looking at the equation (11). For example, one can notice 
that if three sets of measurements are conducted with magnetic induction pointing respec- 
tively in the directions of canonical basis vectors e±, e^, and e^, one can easily reconstruct 
the sum of integrals in the brackets in (11). Further, if one focuses (p(x,t) so that at the 
moment t — it becomes the Dirac 5-function centered at y, i.e. 

<p(x,0) = S(x - y) 

then one immediately obtains the value of V x Jj at the point y (such a focusing is theoreti- 
cally possible as explained in the next section). Thus, by moving the focusing point through 
the object, one can reconstruct the curl of Jj in all of Q. 

Our model also explains the observation reported in [13] that no signal is obtained when 
the acoustic wavepacket is passing through the regions of the constant a(x). In such regions 
current Ji is a potential vector field and, therefore, the integral in (12) vanishes. 

Finally, it becomes clear that an accurate image reconstruction is impossible if monochro- 
matic acoustic waves of only a single frequency k are used for scanning, no matter how well 
they are focused. In this case the spatial component ip of ip(x, t) = if)(x) exp(ikt) is a solution 
of the Helmholtz equation 

Aip + k 2 ip = 0, 

and, within Q it can be approximated by the plane waves in the form exp(?A • x) with 
|A| = k. Let us assume for simplicity that the electrical boundary is removed to infinity. 
Then, measuring M(t) given by equation (12) is equivalent to collecting values of the Fourier 
transform of V x Jj (x) corresponding to the wave vectors A lying on the surface of the sphere 
| A | = k in the Fourier domain. The spatial frequencies of function V x Ji(x) with wave vectors 
that do not lie on this sphere cannot be recovered. 

2 Solving the inverse problem of MAET 

The first step toward the reconstruction of the conductivity is to reconstruct currents Ji(x) 
corresponding to certain choices of I(z). Let us assume that all the measurements are 
repeated three times, with magnetic induction B pointing respectively in the directions of 
canonical basis vectors e±, e 2 , and e 3 . Then, as mentioned above, if ip(x,0) = 5(x — y), one 
readily recovers from the measurements M(0) the curl of the current at y, i.e. V x Ji{y). 
Generating such a velocity potential is possible at least theoretically. For example, if one 
simultaneously propagates plane waves ip\(x,t) = exp(i\ ■ x — i\\\t) with all possible wave 
vectors A, the combined velocity potential at the moment t — will add up to the Dirac delta- 
function S(x). Such an arrangement is unlikely to be suitable for a practical implementation: 
firstly, the sources of sound would have to be removed far from the object to produce a 
good approximation to plane waves within the object. Secondly, the sources would have 
to completely surround the object to irradiate it from all possible directions. Finally, all 
the sources would have to be synchronized. A variation of this approach is to place small 
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point-like sources in the vicinity of the object. In this case, instead of plane waves, spherical 
monochromatic waves or propagating spherical fronts would be generated. These types of 
waves can also be focused into a delta-function (some discussion of such focusing and a 
numerical example can be found in [19]). 



2.1 Synthetic focusing 

However, a more practical approach is to utilize some realistic measuring configuration (e.g. 
one consisting of one or several small sources scanning the boundary sequentially), and 
then to synthesize algorithmically from the realistic data the desired measurements that 
correspond to the delta-like velocity potential. Such a synthetic focusing was first introduced 
in the context of hybrid methods in [18, 19, 24]. It was shown, in applications to AET and to 
the acoustically modulated optical tomography, that such a synthetic focusing is equivalent 
to solving the inverse problem of TAT. The latter problem has been studied extensively, and 
a wide variety of methods is known by now (we will refer the reader to reviews [17, 33] and 
references therein). The same technique can be applied to MAET, as explained below. 



2.1.1 Measuring functionals solve the wave equation 

Let us consider a spherical propagating front originated at the point y. If the initial conditions 
on the pressure p y (x,t) are 

p y (x,0) = 5(x-y), 

!p„M) = o, 

then p y (x,t) can be represented in the whole of M 3 by means of the Kirchhoff formula [32] 

d5(\x-y\-ct) 
ot Atc\x — y\ 

Such a front can be generated by a small transducer placed at y and excited by a delta-like 
electric pulse; such devices are common in ultrasonic imaging. 

Velocity potential <p(x,y,t) corresponding to p y (x,t) then equals 

S(\x-y\ - ct) 

(P{x, V, t) = ——, ; — • 13 

4ti\x — y\ 

The role of variables x and y is clearly interchangeable; ip(x,y,t) is the retarded Green's 
function of the wave equation [32] either in x and t, or in y and t. Moreover, consider the 
following convolution H(y,t) of a finitely supported smooth function h(y) with ip 



H( V ,t) =fm s °' 7 ^'T' ^ 

J An\x-y\ 



Then H(y,t) is the solution of the following initial value problem (IVP) in IR 3 [32]: 

^H(y,t) = A y H(y,t) 



H(y,0) = 0, (14) 

dt J 



f t H(y,0) = h(y). 
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Suppose now that a set of MAET measurements is obtained with propagating wave fronts 
tp(x,y,t) with different centers y (while I(z) and B are kept fixed). By substituting (13) 
into (11) we find that, for each y, the corresponding measuring functional M I<B (y,t) can be 
represented as the sum of two terms: 



M I My,t) = Mf^(y,t) + M^ B (y,t), 



where 



M tjM = B ■ X < Z ) dA ^ (15) 

an 

M^iv, t) = - f S (\*-V\-<*) B . V x Jr(x)dx. (16) 
p J 4n\x — y\ 
n 

It is clear from the above discussion that both terms Mf B g (y,t) and M T T e B (y,t) solve the 
wave equation in M 3 , subject to the initial conditions 

|-Mf° s (a;, 0) = -B- J:{x) x n(x)5 dn (x), (17) 
ot ' p 

^M^ B (x,0) = - p B-VxJj(x), (18) 

Mf^(y,t)=M^ B (y,t) = 0, (19) 

where 5dn(x) is the delta-function supported on dQ. While singular term M s ^ B g {y,t) solves 
the wave equation in the sense of distributions, the regular term Mj° B (y,t) represents a 
classical solution of the wave equation. 



Proposition 1 Suppose conductivity cr(x) and boundary currents I(z) are C°° functions of 
their arguments, and the boundary dQ is infinitely smooth. Then the regular part Mj e B (y,t) 
of the measuring functional is a C°° solution of the wave equation 

^M^(y,t) = A y M^ B (y,t), y G 1R 3 , fe[0,oo), (20) 
satisfying initial conditions (18) and (19). 

Proof. Under the above conditions, potential w i{x) solving the boundary value problem (5), 
(6) is a C°° function in Q due to the classical estimates on the smoothness of solutions of 
elliptic equations with smooth coefficients [12]. Therefore, the right hand side of (18) can be 
extended by zero to a C°° function in IR 3 . Term M] e %{y,t) defined by equation (16) solves 
wave equation (20) (due to the Kirchhoff formula, see [32]) subject to infinitely smooth initial 
conditions (18), (19), and thus it is a C°° function for all y G M 3 , t G [0, oo). ■ 
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2.1.2 Reconstructing the initial conditions 

We would like to reconstruct the right hand side of (18) (and, possibly that of (17)) from 
the measured values of Mj^?/, £). Since c is assumed constant, in the 3D case Mi^{y,t) 
will vanish (due to the Huygens principle) for t > t max = D max / c where -D max is the maximal 
distance between the points of fl and the acoustic sources. We will assume that Mj t B(y,t) 
is measured for all t G [0, t max \. 

The problem of reconstructing -^M I:B (x, 0) from M IjB (y,t) is equivalent to that of re- 
constructing the initial value h(y) from the solution of IVP (14). The latter, more general 
problem, has been studied extensively in the context of TAT (see, e.g. [17, 33] and references 
therein). In the case of TAT, points y describe the location of detectors rather than sources, 
but the rest of the mathematics remains the same. The most studied situation is when 
the detectors are placed on a closed surface £ surrounding the object. If £ is a sphere, a 
variety of inversion techniques is known, including (but not limited to) the explicit inversion 
formulae [11,20,28], series techniques [21,23,27,35], time reversal by means of finite differ- 
ences [2,7, 14], etc. If £ is a surface of a cube, one can use the inversion formula [22], the 
already mentioned time reversal methods, or the fast algorithm developed in [21] for such 
surfaces. 

The choice of the TAT inversion method for application in MAET will depend, in par- 
ticular, on the mutual location of the electric boundary dfl and the acoustic source surface 
S. One can decide to move the electric boundary further away by placing the object in a 
liquid with conductivity equal to 1 and by submerging the acoustic sources into the liquid, 
in which case £ will be inside dfl. Alternatively, one can move the acoustic surface further 
away so that it surrounds the electric boundary. And, finally, one may choose to conduct 
the electrical measurements on the surface of the body, and to place the acoustic sources on 
the same surface, in which case £ will coincide with dfl. 

If dfl lies inside £, all the above mentioned TAT inversion techniques (theoretically) will 
reconstruct both Mj 1 B s (y,0) and Mj° B (y,0). In practice, accurate numerical reconstruction 
of the singular term Mf B g (y, 0) supported on dfl may be difficult to obtain due to the finite 
resolution of any realistic measurement system. Luckily, the support of M] e B (y, 0) (fli in our 
notation) is spatially separated from dfl, so that the contributions from the singular term 
can be eliminated just by setting the reconstructed image to zero outside fli. As explained 
further in the paper, M] c B {y, 0) contains enough information to reconstruct the conductivity. 
On the other hand, M^ B g (y,0) does carry some useful information, which can be recovered 
by a specialized reconstruction algorithm. 

If £ lies inside dfl (but fli still lies inside £), not all TAT reconstruction techniques can 
be applied. Most inversion formulae will produce an incorrect result in the presence of the 
exterior sources (such as the term Mf n B {y, 0) supported on dfl and exterior with respect to 
the region enclosed by £). On the other hand, formula [22], series methods [21] and all the 
time reversal techniques automatically filter out the exterior sources and thus can be used 
to reconstruct Mj e ^(y,0). 

The situation is more complicated if surfaces E and dfl coincide. However, those methods 
that are insensitive to sources located outside dfl are also insensitive to the sources located 
on dfl, particularly to those represented by Mj 1 B g (y,0). These methods can be used to 
reconstruct M] c B (y,Q). In other words, the following proposition holds: 
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Proposition 2 If values of measuring functional M IjB (y,t) are known for all y e X and 

t G [0,-D max /c] 7 i/ie term Mj C B (y,0) can be exactly reconstructed in Q± (by using one of the 
above-mentioned TAT algorithms). Moreover, if the conditions of Proposition 1 are satisfied, 
the reconstruction is exact point-wise. 

Reconstructing the curl In order to reconstruct the curl of the current Ji(x), we need 
to repeat the procedure of finding M^^x, 0) times, with three different orientations of B : 
_B0) = \B\ej, j = 1,2,3. As a result, we find the projections of V x Ji(x) curl on e±, e 2 , e 3 
and thus obtain: 

j=i 

(outside Qi the curl of Ji(x) equals since the conductivity is constant there). 

If, in addition, dQ lies inside £ and Mf^(y,0) has been reconstructed, we obtain the 
term 

Jr(x) x n(x)S ga (x) = ^ £ e~M^ B % (x, 0). (22) 

I I j= i 

2.2 Reconstructing the currents 

The considerations of the previous section show how to recover from the values of the mea- 
suring functionals Mj B u) (y, t) the curl of the current Jj and, in some situations, the surface 
term (22). The next step is to reconstruct current J/ itself. 

2.2.1 General situation 

Let us start with the most general situation and assume that only the curl C = V x Jj 
has been reconstructed. Since current Ji is a purely solenoidal field, there exists a vector 
potential K(x) such that 

Jj(x) = V x K(x) + V(x), (23) 

where K(x) has the form 



n 

and ty(x) is both a solenoidal and potential field. Then there exists harmonic ip(x) such that 

y(x) = V^(ar). (24) 

We know that 

8 

= I- (25) 



J i ■ n \on = ~ w i 



d_ 

On 



on 

Therefore, by combining equations (23) and (25) one obtains 



an 
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and ip(z) now can be recovered, up to an arbitrary additive constant, by solving the Neumann 
problem 

Aif>(x) = 0, x e n 

m ^ , , ;k) ( 26 ) 

dn^^ ) V I I v J 47r(z— j/) 



^(--) = / -(vxJV ;,). :ei«] 



Now Ji(x) is uniquely defined by the formula 

Jj( x ) = V x / —^—dy + V^(x), igSI. (27) 
y 47r(x-y) 
n 

Proposition 3 Under smoothness assumption of Proposition 1 current Ji(x) is given by the 
formula (27), where function ip(x) is the (classical) solution of the Neumann problem (26). 

2.2.2 Other possibilities 

If in addition to the curl V x J/(x), the surface term (equation (22)) has been reconstructed, 
there is no need to solve the Neumann problem. Instead, function ty(x) is given explicitly 
by the following formula [31]: 



My)xn(y) 
4ir{x-y) yy ' 



an 

The final expression for current Ji can now be written as 



Ji(x) = V x x 



C(y) dy+ [ My)xn(y) dA{y) 



An(x — y) J An(x — y) 
n an 



n 



C(y) + J f (y) x n{y )5 m (y) 
4n(x — y) 



W x x J ^ yy> ' '[\ y >" ouyy > dy, 16(1, (28) 



where Q is the closure Q. The term with the delta-function in the numerator of (28) coin- 
cides with the surface term given by equation (22). In order to avoid the direct numerical 
reconstruction of the singular term, one may want to try to modify the utilized TAT recon- 
struction algorithm so as to recover directly the convolution with contained in equation 
(28). The practicality of such an approach requires further investigation. 

Finally, in certain simple domains one can find a way to solve the equation 

C = V x J 7 

for J i in such a way as to explicitly satisfy boundary conditions (25) and thus to avoid the 
need of solving the Neumann problem (26). One such domain is a cube; we present the 
corresponding algorithm in Section 3. 
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2.3 Reconstructing the conductivity 

In order to reconstruct the conductivity we will utilize three currents J^ k \ k = 1,2,3, 
corresponding to three different boundary conditions Ik- As we mentioned before we are 
using three different magnetic inductions j = 1,2,3. As a result, we obtain the values 
of the following measuring functionals Mj B (j)(y, t) 

M IktB u)(y,t) = J I k (z)u U) (z,t)dA(z), i = 1,2,3, A; = 1,2,3, (29) 
an 

where U(j)(z,t) is the electric potential corresponding to the acoustic wave with the velocity 
potential tp(x,z,t) propagating through the body in the presence of constant magnetic field 
B®. Notice that the increase in the number of currents does not require additional 
physical measurements: the same measured boundary values of u^)(z, t) are used to compute 
different measuring functionals by changing the integration weight Ik{z) in equation (4). 

For each of the currents we apply one of the above-mentioned TAT reconstruction 
techniques to compute 0), The knowledge of the latter functions for B^\ j = 

1,2,3, allows us to recover the curls C {k) = V x jW, k = 1,2,3, (equation (21)) and, 
possibly, the surface terms (22). Finally, currents are reconstructed by one of the 
methods described in the previous section. 

At the first sight, finding a(x) from the knowledge of = aVw Ik) k = 1,2,3, is a 
non- linear problem, since the unknown electric potentials Wi k depend on a(x). However, as 
shown below, this problem can be solved explicitly without a linearization or some other 
approximation. Indeed, for any k = 1, 2, 3, the following formula holds: 

= V x — = (v-) x J« + -C^ = ~ (W) x + -C^ 
a \ a J a a 2 a 

so that 

Vlna(x) x J (k \x) = Cf\x), x e Q, k = 1,2,3. 

Now one can try to find V In a at each point in Q by solving the following (in general) 
over-determined system of linear equations: 

Vina x JW = 

Vlnax J( 2 ) =C^ . (30) 
Vlnax J( 3 ) =C^ 

Let us assume first, that currents J^(x), I = 1,2,3 form a basis in M 3 at each point in Vt. 
There are 9 equations in system (30), whose unknowns are the three components of Vina, 
but the rank of the corresponding matrix does not exceed 6. In order to see this, let us 
multiply each equation of (30) by I — 1, 2, 3. (Since the three currents form a basis, 
this is equivalent to a multiplication by a non-singular matrix). We obtain 

f V\na-(JW x J( 2 )) = C«- J( 2 ) 
Vln(7.(J«x J( 2 )) = -C^ • J« 
V\na-(JW x J( 3 )) = CM . JO) 
i Vln(T-(J«x J©) = -C( 3 ) • J« • 
V\na-(JW x J©) = C( 2 ) • J( 3 ) 
Vlna-(j( 2 )xj( 3 )) = -C( 3 )-j( 2 ) 
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In the case of perfect measurements the right hand sides of equations number 2, 4, and 6 
in the above system would coincide with those of equations 1, 3, and 5, respectively, and 
therefore the even-numbered equations could just be dropped from the system. However, in 
the presence of noise it is better to take the average of the equations with identical left sides, 
which is equivalent to finding the least squares solution of this system. We thus obtain: 



Vina - (jW x J( 2 )) 
Vlna-(J« x J( 3 )) 
Vlna-(J( 2 ) x J( 3 )) 



I(C(1) . J(2) _ C (2) . J(l)) 

i(C(i) . j(3) _ C m . j(i)) 

1(C(2) . J(3) _ C (3) . J(2)) 



(31) 



After some simple linear algebra transformations (see Appendix) the solution of (31) can be 
written explicitly as follows: 



Vina 



2J« • (JP) X J(3)) 



(7(2) . J(3) _ (-7(3) . j(2) 

;,/ I _c(i) . jO) + C (3) . j(i) 

C(l) • J(2) _ C( 2 ) • J« 



(32) 



where M is (3 x 3) matrix whose columns are the Cartesian coordinates of the currents 
J«\l = 1,2,3: 



M = J« 



J(2) 



J(3) 



(33) 



Since, by assumption, currents J® form a basis at each point of Q, the denominator in (32) 
never vanishes and, thus, equation (32) can be used to reconstruct V In a in all of fl. Finally, 
we compute the divergence of both sides in (32): 



A In a = -V • 

2 



J« • (J(2) X J(3)) 



(7(2) . J(3) _ (7(3) . j(2) 

,1/ | _C(1) • JO) + C( 3 ) • J« 
C(l) . J(2) _ C (2) . J(l) 



(34) 



and solve the above Poisson equation for In o in Q subject to the Dirichlet boundary condi- 
tions 

Man = 0- ( 35 ) 

The above reconstruction procedure works if currents J^ 2 \ and are linearly 
independent at each point in Q. For an arbitrary conductivity a this cannot be guaranteed. 
There exists a counterexample [8] describing such a conductivity for which a boundary 
condition can be found such that the corresponding current vanishes at a certain point 
within the domain. Clearly, while such a situation can occur, it is unlikely to occur for an 
arbitrary conductivity a, and our method should still be useful in practice. 

Moreover, the condition of the three currents forming a basis at each point in space can 
be relaxed. Below we show that if only one of the currents J^ l \ J^ 2 \ and vanishes (say, 
j(3) — o) at some point and the two other currents are not parallel, the following truncated 
system is still uniquely solvable: 



Vina x jW = 
Vlnax J( 2 ) =C^ 



(36) 
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Indeed, let us multiply via dot product the above equations by and respectively, 
and subtract them. We obtain 

Vina • (J« x J^) = l - (C« • - C< 2 > ■ J«) . (37) 

Now, multiply the first equation in (36) by x . The left hand side will take the form 

(V In a x J«) • ( JW x J^) = V In a • [j« x ( J« x j( 2 ))] 

= V In a • [( J« • J^) J« - ( J« ■ J«) JW] , 

which leads to the equation 

Vina - [(JW • / 2 )) J« - (J« • JW) J( 2 )] = C« • (J« x j( 2 )). (38) 
Similarly, by multiplying the second equation in (36) by (J^ x J^) we obtain 

Vina • • JW) J™ - (JW ■ JW) J«] = ■ (J^ x J«). (39) 

Equations (37), (38) and (39) form a linear system with three equations and three unknowns. 
In other to show that the matrix of this system is non-singular, it is enough to show that 
the vectors given by the bracketed expressions in (38) and (39) are not parallel. The cross- 
product of these terms yields 

[...] X [...] = (J« X J&) [(J« • J( 2 )) 2 - (J« • JW) (J^ . JO) 

The above expression is clearly non-zero if and are not parallel, and therefore the 
system of the three equations (37), (38), (39) is uniquely solvable in this case. 

Theorem 4 Suppose that the conditions of Proposition 1 are satisfied, and that the conduc- 
tivity cr(x) and boundary currents Ik, k = 1,2,3, are such that at each point x £ Q two of 
the three correspondent currents are non-parallel. Then the logarithm of the conduc- 
tivity In a is uniquely determined by the values of the measuring functionals Mj B u)(y,t), 
k = 1, 2, 3, j = 1, 2, 3, y G E and t e [0, D nmx /c]. ' 

Proof. By Propositions 1, 2, 3, and 4, from the values of Mj B a)(y,t) one can reconstruct 
currents J <yk \x), k = 1,2,3, at each point x in f2. Since at each point at least two of the 
three currents (lets call them jW, and J^) are not parallel, system of the three equations 
(37), (38), (39) is uniquely solvable, and V In a can be found at each point in f2. Since it was 
assumed that a is bounded away from zero, V In a is a C°° function in Q. By computing the 
divergence of V In a the problem of finding In a reduces to solving the Poisson problem with 
zero Dirichlet boundary conditions in a smooth domain Q. m 

The condition that out of the three currents two are non-zero and not parallel, is less 
restrictive than the requirement that the three currents are linearly independent. To the best 
of our knowledge, there exists no counterexample showing that the former condition can be 
violated , i.e. that one of the three currents generated by linearly independent boundary 
conditions vanishes and the other two are parallel at some point in space. On the other 
hand, we know of no proof that this cannot happen. 
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3 The case of a rectangular domain 



In the previous section we presented a theoretical scheme for finding the currents J^ k \ 
k — 1, 2, 3, and conductivity a(x) from the MAET measurements in a rather general setting, 
where the electrical domain Q and surface E supporting the acoustic sources are quite ar- 
bitrary. This scheme consists of several steps including the solution of the inverse problem 
of TAT in the domain surrounded by E, solution of the Neumann problem for the Laplace 
equation in Q, and solution of the Poisson equation in Q for lncr. All these problems are 
well-studied and various algorithms for their numerical solution are well known. However, 
in the simplest case when Q is a cube and E coincides with dQ, the reconstruction can be 
obtained by means of an explicit series solution as described below. 



3.1 Fast reconstruction algorithm 

Let us assume that the domain Q is a cube [0,1] x [0,1] x [0,1], and the sound sources 
are located on dQ. (In practice such a measuring configuration will occur if the object is 
placed in a cubic tank filled with conductive liquid, and the sound sources and electrical 
connections are placed on the tank walls). We will use three boundary conditions I k defined 
by the formulae: 

{\, x E dQ, x k = 1 
-\, x G dQ, Xk = , A; = 1,2,3, x = (xi, x 2 , x 3 ). (40) 
0, x G dQ, < Xk < 1 

As before, all the measurements are repeated with three different direction of the magnetic 
field = \B\ej, j = 1,2,3, and the values of functionals M t B (j)(y,t) (see equation 
(29)) are computed from the measurements of the electrical potentials u^\x,t) on dQ, for 



t G 



0,^ 



We start the reconstruction by applying to Mj B (j){y,t) the fast cubic-domain TAT 
algorithm [21] to recover the regular terms M 1 ^ U) (x, 0), j = 1,2,3, k = 1,2,3. (The 
algorithm we chose automatically sets to zero the sources corresponding to the surface terms 
j%M*™f bU) (x, 0) supported on E). This is done for all three directions of so that we 

immediately obtain the curls C (fc) (x) = V x J^ k \x), k — 1, 2, 3, x G dQ. 
Now, since the currents are divergence- free, 

V x C {k \x) = V x V x J (k \x) = -Aj( k \x), 

and we can try to solve the above equation as a set of Poisson problems for the components 
j( k \x) in Q. Below we discuss how to enforce the correct boundary conditions for these 
Poisson problems. 
Let us recall that 

J^(x) = aVw^(x), 

where w^ k \x) is the corresponding electric potential. It is convenient to subtract the linear 
component of potentials, i.e. to introduce potentials w^'°(x) and J^'°(x) defined as follows: 

w {k \x) =w {k) '°(x)+x k , 

J (k \x) = J {k) '°(x) +e k , fc = l,2,3. 



15 



Now each w^'°(x) satisfies zero Neumann conditions on <9f2, and thus can be extended by 
even reflections to a periodic function in M 3 . Since w^'°(x) is a harmonic function near 
dfl, such an extension would be a harmonic function in the neighborhood of dQ and its 
reflections. Therefore, each w^'°(x) can be expanded in the Fourier cosine series in f2, and 
the derivatives of this series yield correct behavior of the so-computed Vw' i: ' ,0 (i) on <9f2, and 
also of currents J^'°(x) (the latter currents coincide with Vw^'°(x) in the vicinity of <9f2). 
Therefore, the components of currents J^'°(x) should be expanded in the Fourier series 
whose basis functions are the corresponding derivatives of the cosine series. In other words, 
the following series yield correct boundary conditions when used as a basis for expanding 
currents k = 1,2,3: 

f oo oo oo 

J[ (x) — E E A m' n sm7r ^ x i COS7rm;r 2 cos7rna;3 

1=1 m=0 n=0 

oo oo oo ^ 

< (x) = E E A,m,n C0S Sm 7r77la; 2 COS 7mX 3 , (41) 

j=0m=ln=0 
(k) oo oo oo 

Jg (^) = E E E A !mn S ^ n 7r ^ 1 COS 7rmx 2 COS 7TOX 3 
k |=0m=0n=l 

where J^°(x) = (j^' ^), J 2 (fe),0 (x), 4 fe) '°(x)) . Now, since V x J^°(x) = V x J«(x): 

AjW'°(x) = -V x C {k) {x), A; = 1,2,3 (42) 

and, if the Poisson problems (42) are solved for each component of using sine/cosine 
Fourier series (41), the correct boundary conditions will be attained. The computation can 
be performed efficiently using the Fast Fourier sine and cosine transforms (FFST and FFCT). 

One can notice that before Poisson problems (42) can be solved, the curl of the C^ k \x) 
needs to be computed. We compute it by expanding C^ k \x) in the Fourier sine series and 
by differentiating the series, again using FFST and FFCT. This computation amounts to 
numerical differentiation of the data that may contain a significant amount of noise. However, 
this does not give rise to instabilities, since this differentiation is immediately followed by 
an inverse Laplacian (describing solution of the Poisson problem), so that the combined 
operator is actually smoothing. 

Finally, once the currents are reconstructed, we form the right hand side of the system 
(34) and solve this Poisson problem for lna(x) by using the Fourier sine series which yields 
the desired boundary conditions (35). Again, FFST and FFCT are utilized here and in the 
computations of the divergence needed to form equation (34). 

One could notice that all the steps of the present algorithm are explicit, and can be 
performed using FFST and FFCT, so that the number of floating point operations (flops) 
required to complete the computations is 0(n 3 In n) for a (n x n x n) computational grid. The 
same is true for the TAT reconstruction technique [21] used on the first step of computations, 
so that the whole method is fast; it has complexity of 0(n 3 lnn) flops. 

3.2 Numerical results 

Since the author does not have at his disposal any real MAET measurements, the work of 
the present reconstruction algorithm (for a cubic domain) will be demonstrated on simulated 
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(a) (b) 

Figure 1: A 3D simulation with a smooth phantom representing lna(x) 

(a) the cross section of the phantom by the plane x 3 = 0.5 

(b) reconstruction from the data with added 50% (in L 2 sense) noise 




Figure 2: The cross section of the reconstructed image (shown in Figure 1) by the line 
xi = 0.25, x% = 0.5. The thick black line represents the phantom, the gray line corresponds 
to the reconstructed image 

data. The most thorough simulation would entail solving equation (2) with the boundary 
conditions (3) for various velocity fields V(x, t) corresponding to the propagating spherical 
fronts originating at different locations on dQ. For a good reconstruction, the number of 
the data points should be comparable with the number of unknowns. We would like to 
reconstruct the conductivity on a fine 3D grid (say, of the size 257 x 257 x 257), which 
implies having several millions of unknowns. Thus, we would need to solve equation (2) 
several million times. This task is too challenging computationally. 

Instead, for a given phantom and for the set of functions 1^ given by equations (40) 
we computed currents by solving equation (5) with boundary conditions (6). Next, 
the wave equation (14) with the initial condition h(y) = ■ V x J^ k \x) was solved for 
j = 1, 2, 3, k = 1, 2, 3; the values of the solution of this equation at points in dQ simulated the 
regular part pMj C ^(x, t) of the measuring functionals pMi^iUit)- Note that for simplicity 
we did not model the singular term Mf^(x,t) given by (15). Theoretically, when the TAT 
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Figure 3: Plot of one of the simulated measurement functional for one of the detector 
positions (see text for details): thick black line represenents accurate values, gray line shows 
the data with added 50% noise 

reconstruction algorithm [21] is applied to such data, since this term is not supported in fl 
it will not contribute to the reconstruction, and so the additional effort to simulate it would 
not be rewarded with additional insight. 




Figure 4: Profile of the reconstructed curl (gray line) compared to the accurate values (thick 
black line). Shown are the values of the third component of along the line x 2 = 0.5, 
x 3 = 0.5 

As the first phantom simulating lna(x) we used the following linear combinations of four 
C 8 radially symmetric functions (p(x — x^) with centers xW lying in the plane x% = : 

4 



/(*) 


= ^anp^x -x [t 
i=i 


} |), 






= (0.25,0.25,0), 


di = 


0.5, 


x^ 


= (0.25,0.75,0), 


a 2 = 


-0.5, 


x ( V 


= (0.75,0.25,0), 


a 3 = 


-0.5, 


s< 4 > 


= (0.75,0.75,0), 


04 = 


0.5, 



where ip(t) is a decreasing non-negative trigonometric polynomial on [0, r ], such that <p(0) = 
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0C-, 



0.25 



0.75 



0.5 



1.0 




1.0 



Figure 5: The second, (almost) piece-wise constant 3D phantom. Shown are cross sections 
by the planes Xj = 0.25, j = 1, 2, 3 




Figure 6: Plot of the values of one of the simulated measurement functionals for one of the 
detector positions (see text for details): thick black line represenents accurate values, gray 
line shows the data with added 100% noise 

1, (p{t) — for t > ro and the first eight derivatives of p vanish at and at r . Radius r was 
equal to 0.34 in this simulation. A gray scale picture of this phantom is shown in Figure 1(a). 
Figure 1(b) demonstrates the cross-section by the plane x 3 = 0.5 of the image reconstructed 
on a 129 x 129 x 129 computational grid from simulated MAET data with added simulated 
noise. The acoustic sources were located at the nodes of 129 x 129 Cartesian grids on each 
of the six faces of cubic domain Q. For each source 223 values of each of the measuring 
functionals were computed, representing 223 different time samples or, equivalently, 223 
different radii of the propagating acoustic front. 

The measurement noise was simulated by adding values of uniformly distributed random 
variable to the data. The so-simulated noise was scaled in such a way that for each time 
series (one source position) the noise intensity in L 2 norm was 50% of the intensity of the 
signal (i.e. of the L 2 norm of the data sequence representing the measuring functional). In 
spite of such high level of noise in the data, the reconstructed image shown in Figure 1(b) 
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contains very little noise. This can also be verified by looking at the plot of the cross section 
of the latter image along the line x 2 = 0.25, presented in Figure 2. 
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(c) 



(d) 



Figure 7: Simulation in 3D with phantom shown in Figure 5; (a) and (c) are the cross 
sections of the phantom by planes x 2 = 0.25 and X\ = 0.25, respectively; (b) and (d) are 
the corresponding cross sections of the reconstruction from the data with added 100% (in 
L 2 sense) noise 

In order to better understand the origins of such unusually low noise sensitivity, we 
plot in Figure 3 a profile of one of the time series, M Il B (3)(y,t) for point y = (1,0.5,0.5). 
The thick black line represents the accurate measurements, the gray line shows the with 
the added noise. In Figure 4 we plot a profile of the reconstructed curl C^(x) (gray line) 
against the correct values (black line). (This plot corresponds to the cross-section of the 
third component of along the line x 2 = 0.5, x 3 = 0.5). The latter figure shows that 
noise is amplified during the first step of the reconstruction (inversion of the spherical mean 
Radon transform). This is to be expected, since the corresponding inverse problem is mildly 
ill-posed, similarly to the inversion of the classical Radon transform. However, on the second 



20 







- 






PI 






- 








1— . 
























- 




























r 




— 




















— ' 




1 — 








- 






— 








_ 










Figure 8: The cross section of the reconstructed image by the line x\ = 0.25, x 3 = 0.25. The 
thick black line represents the phantom, the gray line corresponds to the image reconstructed 
from the data with added 100% (in L 2 sense) noise 

step of the reconstruction, corresponding to solving the problem (34), the noise is significantly 
smoothed out. This is not surprising, since the corresponding operator is a smoothing one. 
As a result, we obtain the low-noise image shown in Figure 1(b). 

The second simulation we report used an (almost) piece- wise constant phantom of lna(:r) 
modeled by a linear combination of several slightly smoothed characteristic functions of balls 
of different radii. The centers of the balls were located on the pair-wise intersections of planes 
X\ = 0.25, x 2 = 0.25, x 3 = 0.25, as shown in Figure 5. The minimum value of lna(:r) in this 
phantom was (dark black color), the maximum value is 1 (white color). The simulated 
MAET data corresponded to the acoustic sources located at the nodes of 257 x 257 Cartesian 
grids on each of the six faces of cubic domain Q. For each source, a time series consisting 
of 447 values for each measuring functional were simulated. In order to model the noise, to 
each of the time series we added a random sequence scaled so that the L 2 norm of the noise 
was equal to that of the signal (i.e. 100% noise was applied). 

In Figure 6 we present the profile of the time series M IlB (3)(y,t) for the point y = 
(1,0.5,0.5). As before, the thick black line represents the accurate measurements, and the 
gray line shows the data with added noise. 

The reconstruction was performed on the grid of size 257 x 257 x 257. The cross sections 
of the reconstructed image by planes x\ = 0.25 and x 2 = 0.25 are shown in the Figure 7(b) 
and (d), next to the corresponding images of the phantom (i.e. parts (a) and (c) of the latter 
figure). The cross section profile of the image shown in part (d), corresponding to the line 
X\ = 0.25, x 3 = 0.25, is plotted in Figure 8. 

As in the first simulation, we obtain a very accurate reconstruction with little noise. This 
is again the result of a smoothing operator applied when the Poisson problem is solved on 
the last step of the algorithm. An additional improvement in the quality of the image comes 
from the rather singular nature of the second phantom. Indeed, while the noise is more 
or less uniformly distributed over the volume of the cubic domain, the signal (the non-zero 
In a) is supported in a rather small fraction of the volume, thus increasing the visual contrast 
between the noise and the signal. 
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Final remarks and conclusions 



Mathematical model 

In Section 1 we presented a mathematical model describing the MAET measurements. In 
general, it agrees with the model used in [13,30]. However, instead of point- wise electrical 
boundary measurements we consider a more general scheme. The advantage of such an 
approach is generality and ease of analysis and numerical modeling. In particular, it contains 
as a partial case the pointwise measurement of electrical potentials (reported in [13, 30]). 

Another novel element in this model is the use of velocity potentials which allow us to 
simplify analysis and obtain a better understanding of the problem at hand. We discussed 
in detail the case of acoustic signal presented by propagating acoustic fronts from small 
sources. However, the same mathematics can be used to model time-harmonic sources. 
Since the problem is linear with respect to the velocity potential, the connection between 
the two problems is through the direct and inverse Fourier transforms of the data in time. 
Finally, plane wave irradiation (considered for example in [30]) is a partial case of irradiation 
by time harmonic sources, when they are located far away from the object. 

General reconstruction scheme 

In Section 2 we presented a general scheme for the solution of the inverse problem of MAET 
obtained under the assumption of propagating spherical acoustic fronts. (As we mentioned 
above, a slight modification of this scheme would allow one to utilize time harmonic sources 
and plane waves instead of the fronts we used). The scheme consists of the following steps: 

1. Apply one of the suitable TAT reconstructions techniques to measuring functionals 
Mj(k) jB u)(y,t), j,k = 1,2,3, to reconstruct the regular terms B(j) (x,t) at t — 
and thus to obtain the curls of . 

2. Compute currents from their curls (this step may require solving the Neumann 
problem for the Laplace equation) 

3. Find Vina at each point in fl using formula (32) or by solving system of equations 
(37), (38), (39). 

4. Find values of A In a by computing the divergence of V In a. 

5. Compute lncr by solving the Poisson problem with the zero Dirichlet boundary condi- 
tions. 

Theoretical properties and numerical methods for all three steps are well known. The 
first step is mildly ill-posed (similar to the inversion of the classical Radon transform), 
the second step is stable, and the third step is described by a smoothing operator. Our 
rather informal discussion suggests that the total reconstruction procedure is stable (it does 
not exhibit even the mild instability present in classical computer tomography), and our 
numerical experiments confirm this assertion. We leave a rigorous proof of this conjecture 
for the future work. 
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Comparison with AEIT 



MAET is similar to AET in that it seeks to overcome the instability of EIT by adding the 
ultrasound component to the electrical measurements. However, MAET has some advan- 
tages: 

1. The arising problem is linear and can be solved explicitly. 

2. The AEIT measurements seem to produce a very weak signal; successful acquisition of 
such signals in a realistic measuring configuration have not been reported so far. The 
signal in MAET is stronger; in fact, first reconstructions from real measurements have 
already been obtained [13]. 

The case of a rectangular domain 

In Section 3 we presented a completely explicit set of formulae that yield a series solution 
of the MAET problem for the case of the cubic domain. It reduces the problem to a set 
of sine and cosine Fourier transforms, and thus, it can be easily implemented using FFTs. 
This, in turn, results in a fast algorithm that requires 0(n 3 \nn) floating point operations 
to complete a reconstructions on a n x n x n Cartesian grid. 

Feasibility of reconstruction using two directions of B 

It is theoretically possible to shorten the potentially long acquisition time by reducing the 
number of different directions of B. If only two orthogonal directions of magnetic field B 
are used, only two components of a curl C = V x J I will be reconstructed on the first step 
of our method (say C\ and C 2 ). However, since divcurl J = 0, 

d d d 

t\ — C 3 = — — — C\ — — — C 2 . 

ox 3 ox i OX2 

Since C vanishes on dQ, the above equation can be integrated in x 3 , and thus C3 can be 
reconstructed from C\ and C 2 . A further study is needed to see how much this procedure 
would affect the stability of the whole method. 

Acknowledgements 

The author gratefully acknowledges support by the NSF through the DMS grant 0908208. 



23 



Appendix 



Consider the following system of linear equations 




(A x B) 
(A x C) 
(B x C) 



Ri 
R 2 



(43) 



where A, B, and C are given linearly independent vectors from R 3 , X G I 3 is the vector of 
unknowns, and Rj, j = 1,2,3, are given numbers. The first equation can be re-written in 
the following form 





Xi 


X 2 


x 3 




X\ 


d\ 


61 




X\ 


h 


Cl 


Ri = 


CLi 


a 2 


a 3 




x 2 


a 2 


b 2 




x 2 


b 2 


c 2 




h 


b 2 


b's 




X3 


a 3 


h 




X3 


h 


c 3 



or 



n = 



detM 



Xi 

x 2 
^3 



61 
b 2 



Cl 

c 2 

C3 



(44) 



where M is a (3 x 3) matrix whose columns are vectors A, B, and C, and r\ = R\j detM. 
Similarly, 



1 



r 2 



detM 



and 



1 



^3 = 



detM 



a 2 

^3 



X\ 

x 2 

X3 

bi 
b 2 
h 



Cl 

c 2 
c 3 

Cl 

c 2 
c 3 



(45) 



(46) 



where r 2 = — i? 2 / detM and r 3 = _R 3 / detM. Formulae (44)- (46) can be viewed as the 
solution of the following system of equations obtained using Cramer's rule: 

ai b\ ci 
a 2 b 2 c 2 
03 h c 3 

Therefore, solution of system (43) is given by the formula 



?"3 ^ 




( Xi 


r 2 


\-\ 


x 2 


ri J 




{ X3 



X = 



detM 

In addition, det M = A ■ (B x C). 



ai bi Ci 
a 2 b 2 c 2 
03 h c 3 
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